Prom force distribution to average coordination number 
in frictional granular matter 
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We study the joint probability distribution of normal and tangential frictional forces 
in jammed granular media, fn(/i,/n), for various friction coefficient ji, especially when 
H — oo. A universal scaling law is found to collapse the data for /i = to oo demon- 
strating a link between force distribution P^{ft, fn) and average coordination number, 
z^. The results determine z^' for a finite friction coefficient, extending the constraints 
counting argument of isostatic granular packing to finite frictional packings. 



Granular matter undergos the jamming transition 
evolving into an amorphous state with a non-zero yield 
stress as the density increases to a point where all par- 
ticles are in contact [H. It has been shown experimen- 
tally and numerically that forces are inhomogeneously 
distributed within a jammed granular system, and fur- 
ther appear to decay exponentially or stretch exponen- 
tially for large values of the force 0, H, 0, [E]- To date, 
there are various theoretical attempts to describe the 
force distribution predicting different behavior. For in- 
stance, lattice models and like Boltzmann-equation ap- 
proaches Q predict an exponential decay. Attempts to 
fit experimental data within the energy ensemble |7| pre- 
dict stretched exponential behavior. But the predictions 
are difficult to justify, since for granular matter energy 
is neither well defined nor conserved due to frictional 
forces. An alternative approach is to use the so-called 
force canonical ensemble with a Boltzmann distribution 
where the boundary stress, not energy, is the conserved 
quantity [1, [^, . It is of interest to reduce the above de- 
fined force ensemble to obtain a single force distribution, 
but methods to accomplish this remain in their infancy 
mainly due to the lack of knowledge on the density of 
states 0. A crude approximation would ignore correla- 
tions between forces and the contact network as well as 
the density of states and would predict an exponential 
decay for the force distribution fg, Q . 

Besides the force distribution and the density of states, 
an additional quantity of interest in this study is the av- 
erage coordination number, z^, of a system at the jam- 
ming transition with interparticle friction coefficient /i. 
Despite the importance of for determining the pack- 
ing stability, there is only one theoretical framework to 
characterize related to the counting argument of the 
isostatic conjecture [ll|. At the isostatic limit, the config- 
uration of contact forces has a unique solution if the con- 
tact network is given, since the number of independent 
forces is identical to the number of balance equations. 
Previous works 0, [H, [H, [13, [H, [H, [13 have shown that 
packings at the jamming transition point are isostatic 
[12j only for two extreme cases, ji — Q and /i = oo, with 
average coordination number zj? — 2d and z'^ = d H- 1 re- 



reaches minimum at = and oo. 

Lacking more definite theoretical approaches to under- 
stand the force distribution, the density of states and 
for a general /i, we perform a numerical study of 
the joint force distribution in frictional granular matter, 
Pfj.{ft,fn)- Here, the forces at the contacts are normal- 
ized by the average forces, in the tangential direction 
/t — and in the normal direction /„ — i^^- We 
show that the key distribution is that of infinite /i, inter- 
preted in terms of the density of states and exponential 
statistics, providing guidance to theoretical attempts un- 
der the statistical framework. We show a universal form 
of the force ratio distribution P^{u), where u is the ratio 
of normal and tangential force, u — jr- valid for all fi, and 
a scaling law is found to collapse all the (u) determin- 
ing zj^ for packings. By using P^{u) we introduce a way 
to calculate the average coordination number for various 
/i based on the Maxwell construction of constraint argu- 
ments. Thus, we extend the isostatic condition from the 
limits of /I = and /i = oo to finite /i, providing the scal- 
ing of z^, an unsolved nonlinear problem. Our results 
provide a connection between two important quantities 
to describe jammed matter: from force distribution to 
coordination number. 

The packings we studied are composed of 10,000 equal 
size spheres interacting with Hertz forces along the con- 
tact direction, F^, and Mindlin forces in the tangential 
direction. Ft , plus the Coulomb condition. Ft < iiF^ [l^ . 
We first generate a gas state without friction at an ini- 
tial volume fraction 0j, then the packing is prepared with 
friction through a slow compression and relaxation pro- 
cess to achieve equilibrium at a given volume fraction and 
coordination number as close as possible to the limiting 
density of the jamming transition. A detailed description 
of the simulation is given in p^j . 

We start by constructing an empirical formula of 
Pooift, fn) based on two numerical results: 

(i) We find that the ratio force distribution [1, 0, [H, 



spectively, where d is the dimension. Recent studies 13 1 



confirm that the indeterminacy of the force ensemble 



P^(u) = K / fnPf^inufn, fri)df„ 



(1) 



14l | at infinite friction is characterized by two power-laws 
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FIG. 1: (a) and the inset are log-log plots of the PDF of u 
respectively in 3D and 2D for various /i; (b) Log- log plot of the 
collapsed -P^(u) for various /i in 3D. The red dash-lines both 
in (a) and (b) are plots oiP^{u) = 

where we use k = 3.80 and 3.43 respectively for 2D and 3D 
from a direct measurement of the simulation; 
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FIG. 2: (a), (b) Contour plotting of Poo(/t,/n) from simula- 
tion results in 2D and 3D respectively; (c), (d) Contour plot- 
ting of the empirical formula Eq. ((2)1 with a = 0.8 in 2D and 
3D respectively; (e), (f) Contour plotting of Po.aift, fn) from 
simulation result in 2D and 3D respectively with fi = 0.3. In 
(a), (b), (e) and (f), we superpose the data from 20 individual 
configurations, each of them contains 10,000 grains. 

with exponents equal to and -3 in 2D, and 1 and -3 in 
3D respectively at u — > and u — s- oo, where k = jp^, 

is an anisotropy parameter. Figure [T^ plots Pfj,{u) for 
various values of /j,, showing that all P^i{u) displays sim- 
ilar behavior having two power-law slopes except for a 
sharp peak at u = fi, due to sliding contacts reaching the 
Coulomb threshold. A correct form of force distribution 
should predict this power-law behavior. 

Notice that previous 2D simulations [1, [l^ have re- 
ported a plateau of Pfj,{u) in the region of [0, /Lt], corre- 
sponding to the first power-law of (u) with the expo- 
nent equal to 0, shown in the inset of Fig. [1^. The second 
power-laws only appears for very large values of /i and 
has not been reported by previous studies. We only show 
Pfj.{u) with > 0.1 for 2D in the inset of Fig. [T^ due 
to the difficulty of preparing disordered 2D monodisperse 
packing at small values of fi. 



(ii) We find that the contour plot of Pooift, fn) follows 
the geometric behavior shown in Fig. [5^ and es- 
pecially in 3D case where Pooift, fn) is symmetric in the 
space of (ft, fn). We will show later on that this symmet- 
ric behavior only occurs at large enough forces in 3D. A 
correct form of the force distribution should predict this 
behavior. 

By fitting our numerical data, we find an empirical 
form of Pooift, fn) for infinite friction, consistent with 
(i) and (ii). We describe it by defining new variables 
/ = Vf'f + fn and 9 = arctan(A), 

Pooif,9)^agie) e~^f, (2) 

where a is a constant, which could be regarded as the 
inverse of the angoricity [1, [^, [l^] . By fitting this distri- 
bution we find 

g{e) = (d- l)(sin6')'*"2cos6', 

which can be regarded as the density of states approxi- 
mately for the force ensemble at ^ = oo. Equation ([2|). 
plotted in Fig. [2t and shows similar pattern to the 
simulation results of Fig. [5^ and We further study 
the contour plot of Pf^ift, fn) a-t fi — 0.3 shown in Fig. 
[2^ and [If. Po.siftifn) displays the same pattern inside 
the Coulomb cone as when jj, — oo. We therefore suggest 
that the study of force distribution for frictional packing 
should focus on packings with fi — co. The density of 
state gi9) describes the probability of the contact forces 
for a single contact to have an angle 9 [we note that 
there is no obvious geometric meaning for 9, which is not 
the angle between the normal and the net contact force: 
9 — arctan(j^) = arctan(K^) ^ arctan(-^)] and indi- 
cates that normal and tangential forces are correlated to 
each other even when there is no Coulomb constraint. 

We define P<i9) as the cumulative probability distri- 
bution of 9 indicating the probability of the contact forces 
for a single contact to have an angle less than 9, such that 
9i0) = We find: 

P^i9) ^ ism9)^-\ (3) 

This simple form P< could lead to a theoretical approach 
to the force distribution since it provides the density of 
states within the statistical mechanics framework . 

Equation Q implies that P^if,9)/gi9) = ae'^f is 
independent of 9. We plot Pooif, 9)/gi9) for various val- 
ues of 9 in Fig. [5^ to further compare with the simulation 
results. We find that all the curves collapse with expo- 
nential tails in the region of / > 1, indicating that the 
empirical form of Eq. ^ captures the main features of 
the force distribution for large forces. Pooif ,9) / gid) ^as 
a peak at / ~ 1 when 9 is small, and exhibits a mono- 
tonic exponential decrease when 9 is close to ^ . This im- 
plies that the probabilities of single forces, Pooif n) and 
Pooift), have different behavior as shown in Fig. [SJd: 
Pooif n) displays a peak at / ~ 1 while Pooift) does not. 
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FIG. 3: (a) Log-linear plot of Pca{f, 9)/g{9) for various value 
of 6 in 3D. All the curves well collapse with a pure exponential 
tail in the region of / > 1. The red dashed-line is function of 
ae-'^f with a = 0.8. (b) Log-linear plot of Poo(/t), Poo{fn) 
and Poifn)- 



This result is consistent with previous experimental stud- 
ies of frictional packings . In Fig. [8)3 we plot the force 
distribution at /i = and compare with fi — 00. We con- 
clude that Poifn) has a stretched exponential tail close to 
Gaussian with a exponent f3 = 1.65 due to local entropy 
maximization [19]. 

By using Eq. ([T]) and Eq. we obtain a ratio force 
distribution Poo(w), shown in Fig. [l] as a red dashed-line 
in good agreement with numerical results. This result 
further confirms that our empirical formula Eq. Q is 
reasonable. 

Further, we show that the ratio distribution is the link 
between the ensemble of forces and the average coordi- 
nation number. We find that Pfj,{u) can be rescaled to 
a single curve (except for the peak at /x), with scaling 
factors equal to A and A^, for the y and x axes, respec- 
tively. We find A = 1 in 2D and A = (z° - z^)/{z';! - z^) 
in 3D, as plotted in Fig. [T]d. In the 2D case, col- 
lapses without scaling, so A = 1. The 3D case is different, 
where we find that A — > 1 when /i — s- 00, so Poo(u) does 
not change after multiplying the scaling factors. The fac- 
tor A diverges at ^ = 0, implying that P^{u) reduces to 
a delta function at ^ = due to the fact that all contact 
forces reach the Coulomb threshold in a pure frictionless 
packing. 

Next, we show that the universal form of Pooiu) deter- 
mines z^ for any /x, hereby extending the isostatic count- 
ing argument from /i = and fi = 00 to finite values of /z. 
From linear counting arguments we know that zj? = 2d 
and = d+1, and we want to interpolate to finite /i and 
obtain z^'. Below, we show that the Maxwell constraint 
arguments based on the number of redundant constraints 
provides the framework to derive zj^. Analysis of the co- 
ordination number of granular packings can be related to 
the Maxwell constraints counting in the rigidity percola- 
tion theory [2(3 |: 



the number of constraints, Nr is the number of redundant 
constraints, and z is the coordination number. At the 
jamming transition, F = 0, resulting in a minimum value 
of z, i.e., z^. Here zdN/2 is equal to the total number of 
unknown force variables for a fixed force network. 

We consider a static packing with both force and 
torque balances, but without any typical constraints of 
translation and rotation. For packings with fi = 00, 
the number of constraint Nc will be equal to the num- 
ber of force balance equation, dN, plus the number of 
torque balance equation, d{d — l)N/2, i.e., iVc(oo) = 

Si + l)N/2. There exists reasonable evidence [l|, 0, 
, [ll 111, Q [M [H, [13 to believe that at the jam- 
ming transition, Nr{oo) = 0, implying a conjecture that 
the Maxwell counting approximation is exact. Therefore, 
z = = d -I- 1. Another important case is at = 0. 
Here the redundant constraints, A^r(O) — d(d — l)iV/2, 
is equal to the number of torque balance equation due 
to the absence of tangential force. Further, we must add 
z{d — l)N/2 extra constraints to A^c(O), corresponding 
to equations of tangential force equal to zero. Ft = 0. 
Therefore, 7Vc(0) = A^c(oo) + z{d-l)N/2 and we obtain 
z = z° = 2d. 

Analyzing intermediate values of /u, is complicated 
since many inequality constrains are created as /Lti^„* — 
Ft* > 0. Calculating z^ becomes a nonlinear problem 
and can be understood as an optimization of an out- 
come based on some set of constraints, i.e., minimizing 
a Hamiltonian of the system, H(Fn,Ft), over a convex 
polyhedron specified by linear and non-negativity con- 
straints. An interesting feature found in previous studies 
is that z^ monotonically decreases from 2d to d -|- 1 with 
increasing fi [1, [H, [H, [13], implying that we can map 
this non-linear problem to a linear one by considering a 
monotonic change in the number of constraints in Eq. 
(|4]) with increasing fi. 

The above analysis suggests to extend the Maxwell 
counting argument Eq. to a system with finite /i 
as: 



F 



■4d 
2 



N - iVe(oo) + [Nr{Q) ' Z^{d - l)iV/2] 77(/i) = 0, 

(5) 

where ri{fi) is an undetermined monotonic function rang- 
ing from 1 to as /i ranges from to cxd. The problem is 
reduced to choosing a functional form for rj^j,). 

To determine 77 (/i), we notice that it should be related 
to the sliding rate of packings, i.e., the ratio of the num- 
ber of the sliding contacts to the number of total contacts 
in a packing, denoted S'(/i). By definition, S{^) is deter- 
mined by P^{u), providing a link between coordination 
number an force distribution: 



F=^N-N, + Nr, 



(4) 



where F is the number of degrees of freedom (or floppy 
modes) satisfying F > 0, TV is the number of grains, Nc is 



S{fi) = 1-1 P^iu)du =1-1 \^Poo{.\u)du (6) 
Jo Jo 

The limiting cases are 5(0) = 1 and S{oo) = 0, and S'(/i) 
has the same monotonic behaviour as ?y(/i). 
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Eq. (7) 

□ (^,=0.40, iV=1024, T=100KPa 
O *;=0.53, iV=1024, cr=100KPa 
_ ^,=0.55, iV=1024, T=100KPa 
V ^,=0 57, iV=1024, cr=100KPa 

*,=0.59, iV=1024, (T=100KPa 
<] (^,=0.61, Af=1024, cT=100KPa 

i^,=0.63, iV=1024, ,T=100KPa 

• ♦,=0.63, JV=10000, o=500KPa 

• ♦ =0.525, iV=10000, D=1500KPa 



10" 



10 



10"' 



10' 



FIG. 4: z'^ versus /i for various initial volume fraction (pi in 
3D. The red solid line is the theory result predicted by Eq. 
(O. In the inset, the red-dash and dash-dot lines are the 
prediction of {z'^ — z'')/2 and (z*^ — z°°)/2, respectively, in 
comparison with simulations. 



While must be a function of S{^), there are many 
choices for the functional relation between both quanti- 
ties. We determine this functional form by fitting the 
simulations. Setting r]{iJ,) = 1 — (1 — S{fi))/X provides 
very good fitting of with simulations both in 2D and 
3D. Substituting into Eq. ([5]), we arrive at a cubic 



equation for in 3D: 



6- zi^ 



9 ? 



-3 = 0. 



(7) 



relations, 
for /i ^ 



It can be shown that Eq. ([7]) predicts two power-law 
" ~ ''c ~ and zj^ — z^ ^ Ai^'^j respectively 
1 /i oo, where a = 2/3 and /? = 2. 
In Fig. m we plot z^ obtained from the cubic Eq. ^ 
and compare with simulation data in 3D. The asymptotic 
predictions of a = 2/3 and /3 = 2 are in good agreement 
with simulation results shown in the inset of Fig. ID It 
is difficult to check the value of /3 due to the difficulty of 
preparing a 3D packing as close as possible to z^ — 4. 
To solve this problem, we prepare larger packings slightly 
above the critical point with a small constant pressure, 
and zj! is replaced by z^ without suffix. This result is 



shown in Fig. [4] with two sets of data for pressure a = 
SOOKpa and a = ISOOKPa. We can see that the power 
law of coordination number is independent of pressure 
even when z'^ is far from the isostatic value. 

When we combine power-law finding of with our 
theoretical work of [17] in 3D, where z^ is linked to the 

zl^ + 



and /z; 



volume fraction (/)^ with a simple formula, (, 
2\/3), then we solve the relation between 
0^ follows the same scaling behavior with /i, (/)° — ^ ~ 
fi", and (j)^ — (jf^ ~ ■ Recent experiments [21| in 
3D investigate the preparation of packings close to the 
random loose packing limit. They find a = 0.51 ± 0.25 
and j3 — 0.89 ±0.16. Their measurement of /3 is far 
away from our prediction which could be due to the same 
reason as us, i.e., the difficulty of preparing packing as 
close as possible to = 4. 

In the 2D case, 'q[y) — S{^) since A = 1, and we have 



2Kfi 



(l + K2/i2)l/2 



/ 



(1 -f K2/i2)l/2 



(8) 



This equation predicts a = 1 and (3 = 2, close to our 
simulation result oi (3 = 1.86. We can not determine the 
value of a from simulation due to the difficulty of prepar- 
ing disordered 2D monodisperse packings for small values 
of /i, and the polydispcrsity of packings may slightly af- 
fect these two indices. Previous simulations [lB| of poly- 
disperse 2D packings have a = 0.7, still close to our pre- 
dictions. 

In summary, we develop a framework to study the con- 
nection between the force distribution and the coordina- 
tion number. Some aspects of this connection remain 
empirical, including the density of states, g{0), and the 
scaling factor A, allowing for the collapse of P^{u) into a 
single curve. Overall, the obtained mathematical forms 
of the density of states, the different force distributions, 
the coordination number and volume fraction, may allow 
for their incorporation into a statistical force ensemble of 
jammed matter [1, [l^]. This may facilitate the solution 
of outstanding open problems such as the prediction of 
the power law scaling of the pressure, the coordination 
number and elastic moduli with the volume fraction near 
the jamming transition [1,]. 
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